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Abstract 

The bandgap and band bowing parameter of semiconductor alloys are calculated with a fast and realistic approach. 
The method is a dielectric scaling approximation that is based on a scissor approximation. It adds an energy shift to 
the bandgap provided by the local density approximation (LDA) of the density functional theory (DFT). The energy 
shift consists of a material-independent constant weighted by the inverse of the high-frequency dielectric constant. 
The salient feature of the approach is the fast calculation of the dielectric constant of alloys via the Green function 
(GF) of the TB-LMTO (tight-binding linear muffin-tin orbitals) in the atomic sphere approximation (ASA). When it is 
applied to highly mismatched semiconductor alloys (HMAs) like Zn Te x Sei_ x , this method provides a band bowing 
parameter that is different from the band bowing parameter calculated with the LDA due to the bowing exhibited also 
by the high-frequency dielectric constant. 
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1. Introduction 

Semiconductor alloying is a very important way of 
bandgap engineering with applications in tailoring the 
materials properties like the electronic stucture or opti- 
cal behaviour. The bandgap of semiconductor alloys can 
be usually described by a quadratic polynomial in con- 
centration x with the quadratic coefficient being called 
the bowing coefficient. Accordingly, the bandgap E g of 
an alloy AB x Ci_ x is described by 

E g (x) = xE g (AB) + (1 - x)E g (AC) -bx(l-x),(l) 

where E g (AB) and E g (AC) are the bandgaps of the bi- 
nary constituents AB and AC. The bandgap bowing co- 
efficient b describes the deviation from linearity and is 
close to a constant for most alloys. 

Highly mismatched semiconductor alloys (HMAs) 
exhibit a large bowing parameter. HMAs are present in 
III-V systems like Ga N x Asi_ x JH], Ga Bi t Asi_ c , and 
Ga Sb v Asi_ x [2]; II- VI compounds like Zn S x Tei_ x 
|H, Zn Te x SeiZM, and Zn O x Tei_ x Jl; and IV-IV al- 
loys like Sn A Gei_ x 0|. Lattice relaxation, localization, 
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charge transfer, and huge bandgap bowing are encoun- 
tered in HMA systems due to large differences in atomic 
sizes and orbital energies, and large lattice mismatch. 

The bandgap bowing of HMAs cannot be described 
well by simple models like virtual crystal approxima- 
tion, at least in the mid-alloy regime of concentrations 
|0]. Other models like the band anti crossing (BAC) 
model [1], the valence-band anticrossing (VBAC) mod- 
els |2|] or combinations of them Jill are able to ex- 
plain the strong bandgap change with concentration x. 
These models are applied in the dilute limit of concen- 
trations, in which one constituent is considered as an 
impurity with composition-dependent coupling to the 
conduction band minimum or valence band maximum. 
However, the BAC model and its extensions completely 
neglect cluster states induced by impurities [9]. Super- 
cell methods overcome this shortcoming but the price 
is large unit cells II 1011 . More efficient approaches are 
based on special quasirandom structures (SQSs) having 
smaller supercells in which the most relevant atom-atom 
correlation functions are similar to those of random al- 
loys ifTTll . 

Standard first-principles band theory is based on the 
Kohn-Sham (KS) scheme lfl2"ll of the density functional 
theory (DFT) 11311 . which is, in principle, exact. The 
DFT is a total-energy method for calculating the elec- 
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tronic ground state. The local density approximation 
(LDA) II 1311 and the generalized gradient approxima- 
tion (GGA) [ 14] provide reasonably accurate results for 
ground-state properties like crystal structure, pressure- 
induced transformations, phonon spectra, magnetism, 
etc. Usually the LDA leads to some overbindings. For 
example, in materials, where narrow bands determined 
by semi-core states (for instance, Zn 3d states in ZnSe 
and ZnTe) contribute to the bonding, the LDA pre- 
dicts these localized semi-core states too high in en- 
ergy. In addition, the LDA and GGA do not describe 
well the long-range part of correlations and have the so- 
called "bandgap problem". Even though the LDA often 
yields good band dispersions for the valence bands, the 
LDA bandgap energy in semiconductors and insulators 
is much smaller than the experimental bandgap. The 
theory that may solve that problem is the Hedin's GW 
approximation (GWA: G = Green function and W = 
screened Coulomb interaction) Jl5, Kjl 12]. The cor- 
rections in the GWA are provided by the self-energy 
that resembles the Hartree-Fock exchange term, which 
is a non-local operator. Fortunately the eigenvectors ob- 
tained with the Hartree-Fock approximation (HFA), the 
LDA, and the GWA are often very close HQ]!. There- 
fore, in many cases one can obtain good results with the 
GWA by solving first the LDA problem, and then taking 
the expectation value of the GWA self-energy over the 
LDA eigenvectors. 

The calculations in the GWA are based on the 
screening of the exchange interaction with the inverse 
frequency-dependent dielectric matrix as opposed to lo- 
cal exchange-correlation potentials in the LDA. The 
computational cost, however, is large due to the eval- 
uation of dielectric matrices, their inversion, etc. Thus, 
the GWA is usually restricted to relatively small sys- 
tems. One way to circumvent these computational hur- 



dles is to use dielectric function models 1 20, 21] and di- 
electric models such as Thomas-Fermi screening in the 
screened-exchange LDA (SX-LDA) methods I22I l23ll . 
Although these schemes are computationally "cheaper" 
than the GWA, only quite recently have SX-LDA meth- 
ods been used to larger systems 11241 12511 . Moreover, 
despite the fact that these methods provide a much bet- 
ter description of the bandgaps than the LDA, it appears 
that it is rather difficult to define the screening constants 
of alloys. 

Generally, in the HMA the bandgap bowing coef- 
ficient is composition dependent, reflecting the strong 
wave-function localization at the band edges (see ||4J] 
and the references therein). In such materials there is a 
strong interaction between the band edge states and lo- 
calized, impurity-like states in alloys. To be explicit, in 



semiconductor alloys like ZnTe. v Sei_ v there is a strong 
interaction between the localized Te level and the va- 
lence band of the host in the Se-rich limit of concentra- 
tions. On the other hand, in the Te-rich limit of con- 
centrations there is a strong coupling between localized 
Se level and the conduction band. Because the bowing 
coefficient is calculated as a variation, i.e., a difference 
that enables some cancellations, the LDA would seem 
appropriate, at least in a zero-order approximation. In 
alloys like ZnTe. t Sei_ v , however, the use of the LDA has 
two major caveats. The first is related to the underesti- 
mation of the bandgap and the height of the conduction 
band. The second comes from the overbinding of the 
cation d states that mix too much with p states of the 
valence band maximum and push these p states upward 
in energy. Thus, in HMA the LDA would rather fail 
to have realistic assessment of not only the bandgap but 
also the band bowing over the whole range of alloy com- 
position x. Below, we present a fast and simple method 
that estimates the bandgap and band bowing coefficient 
with an accuracy approaching that of more demanding 
methods like the GWA and SX-LDA. The concrete ex- 
ample of ZnTe t Sei_ v is analyzed and discussed. 

2. Method 

There is still a simpler method to calculate the semi- 
conductor bandgaps. It is based on the scissor operator 
lfl6ll and is basically a dielectric scaling approximation 
(DSA) for correcting the LDA bandgaps 12611 . Thus for 
a wide variety of materials the main correction to the 
LDA bandgap is JH 



A si — , 
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where a = 9.1 eV is a universal constant and is 
the high-frequency dielectric constant of the material. 
The method is applied not only to usual zinc blende or 
wurtzite types but also to other structures like techno- 
logically interesting high-dielectric materials & Be- 
side the technical arguments presented in the original 
paper [26] there are other arguments that ensure the 
success of this method. First, there is the perturbative 
argument of the GWA with respect the LDA 1181 : the 
GWA corrections are basically the expectation value of 
the GW self-energy over the LDA eigenvectors. Sec- 
ond, the difference between quasiparticle energies and 
the LDA eigenvalues is essentially due to the non-local 
nature of the effective potential 12811 . which can be ex- 
pressed with a scissor-like operator N26I1 . There is also 
a more physical picture provided by Harrison [29]. The 
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exact non-local exchange operator of the Hartree-Fock 
theory has a long-range part which is unscreened, and 
hence overestimated. In the LDA the exchange opera- 
tor is replaced by a local exchange-correlation operator 
which neglects the long-range part of the exchange al- 
together. In semiconductors the long-range part of the 
exchange is screened by the high-frequency dielectric 
constant. The LDA uses a fixed one-electron potential 
for all bands including the unoccupied bands. Thus the 
LDA neglects the electrostatic energy U needed for the 
separation of the electron and hole in the excitation that 
provides the bandgap. Harrison noted that the excita- 
tion process is accompanied by a dielectric relaxation 
around the electron and the hole such that the correction 
to the bandgap is a rigid shift for the entire conduction 
band of order U/Soo- Finally, a quite similar approach 
has been reported recently. It starts from the optimized 
effective potential method with a free parameter, whose 
value is established by the separation between the short- 
ranged and long-ranged exchange in the screened hybrid 
functional [30]. 

In the present work we calculate the bandgap and 
band bowing coefficient in semiconductor alloys beyond 
the LDA by calculating the high-frequency dielectric 
constant. The equilibrium positions of atoms in the SQS 
supercells representing the alloy are determined with 
the Vienna ab initio simulation package (VASP) 13 111 
with the projector augmented-wave method (PAW) 113211 . 
Then the electronic structure is calculated with the tight- 
binding (TB) linear muffin-tin orbital method (LMTO) 
in its atomic sphere approximation (ASA) lf33ll34 , 35 1 
in the LDA. The TB-LMTO-ASA method is one of the 
most efficient, versatile, and also one of the most phys- 
ically transparent methods. Other methods use a big- 
ger number of basis functions and they are usually more 
difficult to interpret physically. Even though it uses the 
shape approximation (the Wigner-Seitz (WS) cells are 
replaced by overlapping WS spheres), the TB-LMTO- 
ASA method is quite successful in the description of 
closely packed solids, whereas empty spheres are used 
in the interstitial region to close pack open structures 
03611 . As a rule, the shape approximation does quite 
well for closely packed systems, in particular for en- 
ergy bands, but not that well for total energy calcula- 
tions. Finally the static polarization x° an d the high- 
frequency dielectric constant B m are calculated with the 
TB-LMTO-ASA. The procedure is laid out below. 

Let us then consider the static screening of a test 
charge in the random phase approximation II 1611 . Con- 
sider a lattice of points (spheres), with the electron den- 
sity in equilibrium. We wish to calculate the response 
(screening) in the electron charge Sqj at site j, induced 



by the addition of a small external potential 6Vj at site 
i. The electron charge 6qj is related to SVj by the non- 
interacting response function^: 



(3) 



X®j can be calculated in the first-order perturbation the- 
ory from the eigenvectors of a non-interacting one- 
electron Schrodinger equation with the Adler- Wiser for- 
malism [[37, 38] or directly via G°, the Green function 
(GF) calculated from the one-electron Hamiltonian of 
the LDA. By linearizing the Dyson equation for the per- 
turbed GF, i. e. 5G = G°6VG ~ G°6VG Q , one obtains 
an explicit representation of x°- in terms of G° 



f i = l -\mY J §dz6G ii {Kz) 
= ±X$dzGl{k,z)G Q /k,z)dV° 



(4) 



In Eq. © the k summation is the summation in the Bril- 
louin zone and the energy-integration contour encloses 
the occupied states. A large part of the computer re- 
sources is allocated to the calculation of^- . The calcu- 
lation of x° by the Adler- Wiser formula would require 
the utilization of large computer resources because it 
needs all the occupied orbitals at once. In contrast, Eq. 
<j4j shows clearly that^ can be calculated piecewisely, 
and hence efficiently once the GF is obtained. Equation 
© is the TB-LMTO-ASA counterpart of the static po- 
larization function defined in real space fl6Tl . Similarly, 
the dielectric function is 



e = 1 - Vx°, 



(5) 



where V is the unscreened Coulomb matrix. The macro- 
scopic dielectric constant is calculated as [39] 



l -Y Jim l/[e w (k)]" 
tf ~ k^o 



N„ 



(6) 



In equation (0 N is the total number of orbitals. The 
indices /3 - (R, I, m) denote the lattice sites R and the 
angular-momentum quantum numbers / and m. Never- 
theless equation (|6]l neglects "local field" effects and it 
represents the upper bound of e M when "local field" ef- 
fects are considered l39ll . 

The GF can be directly calculated within the LDA of 
the TB-LMTO-ASA by using the potential parameters, 
C, A, and y JH [HQ 



G(Z) = 



y/P a (z) 



(7) 
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Figure 1: Macroscopic dielectric constant estimated with the TB- 
LMTO-ASA. There is an estimated bowing of about -0.89 



where a is the representation, determined by the screen- 
ing constants cvri and 



■ Cri 



Ar/ + (jri - a R i) (z - Cri) 



(8) 



is the potential function. P a and P a are, respectively, 
the first and the second derivatives of the potential func- 
tion with respect to energy. S a = S" „, is the ma- 

lv/m,lv I'm' 

trix of structure constants in the a representation. In the 
LMTO-AS A approach, the effective potential is entirely 
described by the potential parameters: the band center 
C, bandwidth A, and the distortion parameter y. The 
potential function provides the boundary conditions for 
the radial Schrodinger equation inside each WS sphere 
1 33 , 34 , 4(J • The fact that both the charge density and 
the effective potential are described by a few atomic 
parameters makes the calculation of the response func- 
tions a very simple task. 



3. Results and Discussions 

For illustration we show the calculations of the 
bandgap corrections in the ZnTe. t Sei_ r alloys. The ran- 
dom zinc blende alloys were modeled as supercells of 
54 atoms (27 cations) with special quasirandom struc- 
tures (SQSs) JH. The supercells have the lattice vec- 
tors (3ai, 3a2, 3a3) with (ai, 82,83) the unit cell vectors 
of the binary constituents. The SQS structures are con- 
structed such that the most relevant atom-atom correla- 
tion functions of the SQS structure are similar to those 
of random alloys flTll . We assumed that the alloy lattice 
constants are determined by the weighted average of the 
lattice constants of the constituents (Vegard's law) 14 lfl . 
In the ASA, the spherical average of the charge in each 



sphere makes the TB-LMTO approach less reliable and 
more cumbersome in calculating the equilibrium posi- 
tions of the atoms in the SQS cells. For that reason the 
atomic positions inside the SQS cells were relaxed us- 
ing VASP with the LDA exchange-correlation potential; 
at relaxed positions, the quantum-mechanical forces on 
each atom were less than 0.03 eV/A. As we have al- 
ready mentioned the atomic sphere approximation as- 
sumes that the whole space is filled with (overlapping) 
spheres and that the volume of the interstitial region 
vanishes. This is achieved by adding empty spheres 
l36ll . The empty spheres must satisfy several criteria: 
the sum of the volume of all spheres should equal the 
volume of all space, the average overlap between the 
spheres should be minimal, and the radii of the spheres 
should be chosen such that the spheres overlap in the 
regions close to the local maxima of the electrostatic 
potential. The positions of the empty spheres are well 
known in the high-symmetry systems like zinc blende 
structures. In relaxed SQS cells the positions and the 
radii of the empty spheres were chosen to satisfy the 
above criteria starting from the initial positions of the 
non-relaxed (higher-symmetry) SQS cells. The Bril- 
louin zone integration was performed on a 4x4x4 mesh 
by a Monkhorst-Pack £-point sampling scheme [42]. 

In Fig. [TJ we plot the calculated dielectric constant 
too. The calculated values of for the binary con- 
stituents are 10 - 20% larger than the experimental val- 
ues and 5-15% larger than the time-dependent density 
functional theory (TD-DFT) calculations 114311 . There 
is a bowing shown also by the macroscopic dielectric 
constant e ro , which is b as -0.89. Equation (f2]i and 
the bowing exhibited by e M show us that the DFT-LDA 
alone cannot reproduce the bowing of the bandgap in 
ZnTe A Sei_ A . 

Fig. [2] shows the bandgaps calculated with GGA- 
VASP and LDA- VASP, with the TB-LMTO-ASA in the 
LDA, and with the DS A defined by ©. The TB-LMTO- 
ASA results were considered with combined correc- 
tions |0] ■ The bandgaps calculated with GGA-VASP 
are larger than the results of LDA- VASP over entire 
range of concentrations. The bandgaps of the binary 
constituents calculated with LDA- VASP are both 1.06 
eV. The LDA ASA gives a bandgap value of 1.07 eV 
for ZnSe and 1.05 eV for ZnTe. The ASA values are 
pretty close to the LMTO full potential values which 
are 1.05 eV and 1.03 eV, respectively [17]. Over the 
whole range of Te concentration x, the LDA-ASA re- 
sults are in quite good agreement with the results ob- 
tained with PAW- VASP, which, in principle, is a full po- 
tential method with no shape approximation. Fig. |2]also 
shows the bandgaps calculated with the DSA, which ex- 
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Figure 2: The bandgaps of ZnTe r Sei_ r calculated with the LDA, the 
GGA, and Eq. |2) at various compositions x . The squares and the 
diamonds indicate the GGA-VASP and the LDA-VASP results, re- 
spectively; the triangles and the stars represent the values of the LDA- 
TB-LMTO-ASA and of the DSA, respectively. The continuous lines 
are the results of the fitting. No spin-orbit coupling was considered. 



hibits much larger bandgaps than both the LDA and the 
GGA. The bandgaps of ZnSe and ZnTe calculated with 
the DSA are 2.38 eV and 2.07 eV, respectively, being 
comparable with the results of the GWA, which, in gen- 
eral, underestimates the bandgaps because it generates 
too much screening due to overestimation of x° 017I1 . 

As we expected, the bandgaps show a bowing, i.e., 
the bandgap deviates from linear interpolation with re- 
spect to composition x. The bowing parameter b is 1 .64 
eV for the VASP-GGA and 1.57 eV for the VASP-LDA. 
The LDA of the TB-LMTO-AS A produces a bowing of 
1 .50 e V, while the DSA estimates the bowing parame- 
ter to 1.71 eV. The bowing parameters predicted with 
methods that go beyond the LDA are larger. The fact 
that the LDA does not reproduce the entire bowing of 
the bandgap is rooted not only in the LDA's inability 
to predict the bandgap but also in the overbinding pro- 
duced by the LDA. For ZnTe, ZnSe, and ZnTe A Sei_ x the 
LDA overbinding pushes the cation d levels closer to the 
top of the valence band. Thus the d levels will mix too 
much with the p states at the maximum of the valence 
band; therefore additional errors in estimating the band 
bowing occur. On the other hand it is known that the 
GWA moves these cation d states in the right direction 
downwards 

Old measurements of the bandgap bowing parameters 
b with and without spin-orbit interaction are presented 
in 14411 . The authors estimated that the bandgap bow- 
ing is b — 1 .266 e V. They also estimated the spin-orbit 
splitting, which shows a concave bowing. Thus their 
estimated bandgap bowing without spin-orbit coupling 
is b — 1.07 eV. Other measurements estimate a bow- 
ing parameter of 1.62 eV (at 5 K) and 1.51 eV (at 300 



K) 14511 . More recent experiments set the bowing of the 
bandgap between 1.25 eV and 1.5 eV at 300 K [4, 46] 
and the bowing of the spin-orbit splitting at -0.33 eV 
J3]. Other measurements at low temperature, namely at 
13 K, indicate a bowing parameter around 2.21 eV 14711 . 

Several general features can be outlined from all ex- 
perimental data and our calculations. More recent ex- 
periments reveal larger band bowing coefficients than 
the older ones. The bowing is larger at low tempera- 
tures. There is also a bowing of the spin-orbit splitting 
but it is negative. The bandgap bowing calculated with 
spin-orbit coupling is larger than the bandgap bowing 
without spin-orbit coupling by one third of the bowing 
parameter modulus of the spin-orbit splitting. Overall, 
our results for bandgap bowing are in rather good agree- 
ment with recent experimental results. 

4. Summary and Conclusions 

To summarize, we have presented a computational 
scheme that realistically estimates the bandgaps and 
band bowing parameters with applicability to semicon- 
ductors and their alloys. The scheme uses the static 
polarization function^ to calculate the high-frequency 
dielectric constant e M in the GF-TB-LMTO formalism. 
The corrections to the LDA bandgaps are made with 
a dielectric scaling approximation (DSA), which esti- 
mates the correction to the LDA bandgap as a universal 
constant divided by the high-frequency dielectric con- 
stant. 

The approach inherits its speed from the fact that it 
works piecewisely to calculate^ via GF formalism and 
its efficiency from the TB-LMTO (the use of the muffin- 
tin orbitals and a minimal basis). 

The limitations of the method were already discussed 
in j26ll . The narrow -bandgap semiconductors are par- 
ticularly difficult to describe with the DSA because the 
other corrections of the GWA to the the LDA bandgap 
become as important as the DSA. Also, often in these 
semiconductors there might appear a wrong ordering of 
the LDA levels around the bandgap with poor results in 
GWA too d. 

We applied the scheme to the ternary alloy 
ZnTe A Sei_ A by calculating the bandgap and band bow- 
ing in a supercell with 54 atoms that mimics the al- 
loying of ZnTe A Sei_ A . The ternary alloy ZnTe A Sei_ A 
exhibits large bandgap bowing. The bandgap bowing 
computed with the DSA is in good agreement with the 
experimental data. It was also shown that the bowing of 
the bandgap in ZnTe A Sei_ A cannot be explained by the 
LDA alone due to inherent bowing shown also by the 
dielectric constant. 
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